# This script produces Figure 3 # 
rm(list=ls())
options(digits=4)
pkg <- c("PanelMatch", "ggplot2", "plm", 
         "matrixStats",
         "foreign", 
         "stringi", "stringdist",
         "DataCombine","lmtest", "multiwayvcov",
         "data.table",
         "dplyr", 
         "stargazer",
         "clubSandwich", "mediation",
         "doBy")
lapply(pkg, require, character.only = TRUE)

# load datasets
proprietary <- read.csv("proprietary.csv")
load("non_p.RData")

# merge 
d_sub1 <- merge(d_sub1, proprietary, by = c("county_id", "time"), all.x = T)


matches_PS_L4 <- PanelMatch(covs.formula = ~ 
                              tour + log_transaction_area_l1 + 
                              log_transaction_area_l2 + log_transaction_area_l3 + 
                              log_transaction_area_l4 + 
  
                              log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                              log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                              log_expd_l1 + log_middle_school_l1 + 
                                 
                              msec2currentsec_l1 + 
                              msec2currentgvn_l1 + 
                              mayor2currentsec_l1 + 
                              
                              mayor2currentgvn_l1 + 
                              msec2currentsec2_l1_missing, 
                            match.missing = F,
                            lag = 4, 
                            unit.id = "county_id",
                            time.id = "time",
                            treatment = "tour",
                            refinement.method = "mahalanobis",
                            size.match = 10,
                            outcome.var = "log_transaction_area",
                            #weighting = FALSE,
                            listwise.delete = TRUE,
                            qoi = "att",
                            data = d_sub1)    



matches_PS_L6 <- PanelMatch(covs.formula = ~ 
                              tour + log_transaction_area_l1 + 
                              log_transaction_area_l2 + log_transaction_area_l3 + 
                              log_transaction_area_l4 + 
                              
                              log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                              log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                              log_expd_l1 + log_middle_school_l1 + 
                              
                              msec2currentsec_l1 + 
                              msec2currentgvn_l1 + 
                              mayor2currentsec_l1 + 
                              
                              mayor2currentgvn_l1 + 
                              msec2currentsec2_l1_missing, 
                            match.missing = F,
                            lag = 6, 
                            unit.id = "county_id",
                            time.id = "time",
                            treatment = "tour",
                            refinement.method = "CBPS.match",
                            size.match = 10,
                            outcome.var = "log_transaction_area",
                            listwise.delete = TRUE,
                            qoi = "att",
                            data = d_sub1)

L4_cbps_att_sizes <- summary(matches_PS_L4$att)$overview$matched.set.size
L6_cbps_att_sizes <- summary(matches_PS_L6$att)$overview$matched.set.size

pdf("output/Figure3.pdf")
hist(L4_cbps_att_sizes[L4_cbps_att_sizes>0], 
     col = "grey", border = NA,
     freq = TRUE,
     xlab = "Size of Matched Set",
     main = "Inspection",
     ylab = ""
     )


white.hist <- hist(L6_cbps_att_sizes[L6_cbps_att_sizes>0], 
                   col = rgb(1,1,1,0),     
                   
                   plot = FALSE, lty = 1)
lines(white.hist$breaks, c(0, white.hist$counts),
      type = 'S', lty = 1)

legend(x = 0, y = 400,  
       legend = c("4 Month Lags", "6 Month Lags"),  
       fill = c("darkgrey", "white"),
       
       xjust = 0,
       
       pt.cex = 1,
       bty = "n", ncol = 1, cex = 1, bg = "white")


mtext("Frequency", side = 2, line = 2.2, outer=FALSE,
      cex = 1)
dev.off()

